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One of the major assumptions in the search for gravitational wave signatures from massive and 
supermassive black hole binaries with LISA, is that these systems will have circularized before en- 
tering the LISA bandwidth. Current astrophysical simulations now suggest that systems could have 
a non-negligible eccentricity in the LISA band, and an important level of eccentricity in the Pulsar 
J.— ^ ' Timing regime. In this work, we use a set of source catalogues from astrophysically motivated mod- 

^vj i els of massive black hole binary formation and assume a one year LISA mission lifetime. Depending 

■ on the model in question, the initial eccentricities in the final year of the inspiral can be as high as 

^>' 0.6 for high mass seeds and 0.8 for low mass seeds. We show that restricted post-Newtonian circular 

templates are extremely inefficient in recovering eccentric binaries, with median optimal signal to 
noise ratio recoveries of approximately 10% for all models considered. This coupled with extremely 
large errors in parameter recovery from individual Markov chain Monte Carlo's demonstrate quite 
OO ' clearly that even to search for binaries with initial eccentricities as low as 10~'*, we will require 

^N , eccentric templates for LISA data analysis. 

O^, I. INTRODUCTION. 

U . 

^W. The inspiral and merger of massive and supermassive black hole binaries ((S)MBHBs) will be one of the brightest 

sources of gravitational waves (GW) for the future space-borne GW detector LISA j|. There is currently a lot of 
effort within the community in the development of search algorithms for these sources. Initially, the majority of 
t> ' the development was conducted using non-spinning restricted post- Newtonian (PN) waveforms @t3i i-^- ^^^ phase 

^^ , was constructed at a high PN order, but the amplitude was kept at the dominant Newtonian order. In fact, these 

2^ ■ circular templates formed the basis for the initial Mock LISA data challenges (MLDCs) '1&\ . More recently, people 
have started to investigate more complex circular waveforms such as those with higher harmonic corrections |lll4l5l | 
and spin [13, [l6l - [l8j . In both cases it was comprehensively demonstrated that the added complexity helped to break 

1^ I' parameter correlations and in some cases improved the estimation of parameters. While these templates are definitely 
a step in the right direction towards using more realistic waveforms, they are still based in the circular approximation. 
Until now, very little attention has been paid to the role of eccentric templates, as it had always been assumed 
that a binary system would circularize before entering the interesting frequency band |19l - [2l| . These works showed 
that given an initial eccentricity bq and semi- major axis ao, the instantaneous eccentricity and semi- major axis were 
related by e/eo ~ (a/ao)^^^^^. This inferred that as the semi-major axis shrunk by a factor of 2, the eccentricity 
5J] ■ shrank by a factor of almost 3. Based on this assumption, it was shown in Ref |21|] that circular templates would be 
Ci ^ quite efficient in resolving binaries with a residual eccentricities. However, this study made a number of assumptions 
that proved critical to this result. The first was that the study assumed typical LIGO sources, and was focused on 
systems that it was assumed would (almost) circularize before reaching the lower frequency cutoff for initial LIGO 
at 40 Hz. Secondly, and perhaps most importantly, there were no PN corrections added to the waveform. The study 
used the standard Peters and Matthews equations [I^, [23 to evolve the semi-latus rectum and eccentricity, and not 
the more recent PN equations for radiative dynamics. Thirdly, the study was made for optimally orientated systems 
which allowed the study to neglect one of the GW polarizations. 

More recently, two separate studies in particular have shown the previous results to be erroneous |22l . |23| . These 
particular studies included PN corrections in both the conservative dynamics and radiation reaction. The overwhelm- 
ing conclusion of both of these studies was that circular binaries were incapable of matching eccentric systems once 
the initial eccentricity was greater than 0.1. We should point out, that these studies were again based on a LIGO 
study where the waveform durations are very short (on the timescale of seconds). However it was highlighted in 
Ref 23] that longer duration signals, such as neutron star - neutron star binaries had lower matches with eccentric 
binaries. This was due to the fact that there is more of a possibility of de-phasing between waveforms if the signal is 
of long enough duration. It is for this very reason that it is important to conduct this investigation for LISA, where 
signals have durations of months to years. 

Furthermore, although the circular orbit ansatz has been widely used by the GW community, the reliability of such 
an assumption is now further questionable from an astrophysical point of view. Since the evolution of MBHBs was 
firstly sketched by Begelman Blandford & Rees [2J|, hundreds of studies have been dedicated to the subject. After 
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the two MBHs, driven by dynamical friction [2^ [2g|, reach the center of the merged system and pair together, the 
new-formed binary needs to transfer energy and angular momentum to the surrounding ambient in order to coalesce. 
Typical pairing lengthscales of LISA binaries (binary masses in the range 10^ — 10^ Mq) are of the order of 0.1-1 pc; 
however GW emission is efficient in driving binaries to the final coalescence only at mpc scales. This pc-to-mpc gap in 
the MBHB evolution goes under the name of 'last parsec problem' [23]. In stellar environments, the MBHB evolution 
proceeds via super-elastic scattering of surrounding stars intersecting the binary orbit [slingshot mechanism, [28[ , and 
the fate of the system depends on the supply of stars available for such interaction [29r;32j. On the other hand, if 
the system is gas rich, torques exerted by a massive circumbinary disk have been proven efficient in shrinking the 
binary down to ~ 0.1 pc [3J,[3J], which is the current resolution limit of dedicated smoothed particle hydrodynamical 
simulations. However, whether viscous angular momentum extraction is efficient all the way down to coalescence is 
questionable [35| . Although there are still open questions about the effectiveness of these dynamical processes [36| , 
both stellar and gas based shrinking mechanisms have proven to be efficient in increasing the binary eccentricity 
|37H4a |. Whether a MBHB is going to be circular or eccentric for GW detection purposes, depends on the amount 
of 'residual eccentricity' left by GW-driven circularization by the time the system enters the LISA window. Using an 
hybrid model to couple scattering of bound and unbound stars and GW emission, Sesana [47| carried a systematic 
study of the residual eccentricity of MBHBs evolving in stellar environments. Assuming standard MBH formation and 
evolution scenarios [4^, |40| he demonstrated that LISA MBHBs are likely to show a significant amount of eccentricity, 
calling for the development of trustworthy eccentric templates for MBHB detection and parameter estimation. 

In this paper we test the capability of restricted PN circular templates in detecting GWs from eccentric MBHBs 
for LISA. As well as detection capabilities, we are also interested in any parameter mismatch as a result of the 
inclusion of eccentricity. Furthermore, there has also been incredible progress made in the field of numerical relativity 
(NR) |50l452l | in the last number of years, providing the GW community with catalogues of merger waveforms for 
both non-spinning and spinning black hole binaries. However, a priority in the field has been to find ways of reducing 
the residual eccentricity in the NR waveforms [53, [SJl . It will be interesting to see what our results mean for both of 
these endeavors. 

For this study, we use catalogs of eccentric MBHBs evolved according to the scheme presented in [43]. This 
provided us with six individual catalogues as our starting point. For each model, the binary reaches a separation 
where dynamical friction is no longer efficient in shrinking the binary orbit. As the dynamics leading to this point are 
highly environment dependent, the eccentricity at this point has to be added by hand at the beginning of simulations. 
Thus, the catalogues were generated for both high and low mass seed black hole binaries with initial eccentricities of 
Co = 0,0.3 and 0.6. 

The paper is organized as follows. In Section |TT] we present the response of the LISA detector in the low frequency 
approximation. Section IIIII contains a description of the waveforms for both circular and eccentric binary systems. 
We then describe the astrophysical model used to produce the source catalogues in Section IIVI before describing the 
setup of the Monte Carlo simulation in Section |Vl In Section IVll we present the results of the Monte Carlo runs. 

Throughout the paper we use the units G = c = 1 . 

II. THE LISA DETECTOR RESPONSE IN THE LOW FREQUENCY APPROXIMATION. 

The response h{t) of the A and E LISA channels to an incoming GW with polarizations /i+ and hy \n the low 
frequency approximation [55| is given by the combination 

hA.E{t) = h+{m)FiEit) + h. {m)FiEit), (1) 

where the phase shifted time parameter ^(i) is defined by 

^{t) =t - R(s sin d cos {a{t) -(/)). (2) 

Here R^ denotes 1 AU, {9, (f) are the sky location of the system and a{t) = 2Tr fmt + n, where /„ is the LISA 
modulation frequency and k is the longitudinal offset of the detector array. The polarizations of the GW h+^y{t) will 
be defined at later stage. The beam pattern functions are defined by 

FaA^ = \ [cos(2V')7^+ ^(i; 0, 0, A) - sin(2^)I?^_^(t; 0, 0, A)] , (3) 

F1,e{^) - \ [sin(2^)^XB(^; ^' '^' ^) + <^os{2^)DXE(t- 9, 0, A)] . (4) 



The quantity ip is the polarization angle of the wave. Formally, if L is the direction of the binary's orbital angular 
momentum and fi is the direction from the observer to the source (such that the GWs propagate in the — n direction), 
then Tp fixes the orientation of the component of L perpendicular to fi. The time dependent quantities D^^x{t) are 
given in the LFA by [Sgj 
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6 sin(6') ( cos(3a(f) - 2A - 0) + 3 cos(a(i) - 2A + . 



where A = and tt/4 give the orientation of the two detectors, thus defining the A and E channels. To ensure that 
we are working in the LFA domain, we limit the waveforms to a GW frequency of fgwiGm) or 5 mHz, whichever is 
lower. 

III. THE GRAVITATIONAL WAVEFORM FROM CIRCULAR AND ECCENTRIC BINARIES. 

For this study, we use circular and eccentric non-spinning waveforms. For the circular templates, we use a waveform 
with 2 PN corrections to both the orbital phase and angular frequency evolution. For the eccentric waveforms we 
assume 2 PN corrections to both the conservative and adiabatic dynamics of the system. To make the comparison with 
the eccentric waveforms, we are using the standard restricted PN circular templates that have been used extensively 
in other studies. The reason for this is, we would like to study the fidelity of the non-spinning circular templates that 
have been adopted by the community in capturing the eccentric systems, rather than using the zero eccentricity limit 
of the eccentric waveforms. We will therefore describe each waveform in turn. 

A. The Circular Binary Waveform 



The GW polarizations for a non-spinning circular binary are given by 
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Here vn = m\ -I- m-i is the total mass of the binary, t] ~ mim2/rn? is the reduced mass ratio and D^ is the luminosity 
distance, which is related to the redshift z of the source by 
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where we use the WMAP values of {VLR,nM,^A) = (4.9 x 10-^ 0.27, 0.73) and a Bubble's constant of Hq=11 
km/s/Mpc [57| . The inclination of the orbit of the binary system is formally defined as cos t = L • n. The quantity ipc 
is the orbital phase constant at coalescence and the invariant PN velocity parameter is defined by x{t) = {muj{t)) ' , 
where 
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is the 2 PN order orbital angular frequency for a circular orbit formally defined as w = d^/dt, and $ is the orbital 
phase which is given by 



m = i 
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We should note that the gravitational wave phase is defined by ^gw = 2$. The time dependent quantity Q{t; tc) is 
related to the time to coalescence of the wave, tc, by 



e(t;i,) = -^(i,-t). 



(12) 



As the PN waveforms are known to break down in certain cases before we reach r = 6m, we use a taper function to 
smoothly truncate the waveform. 



B. The Eccentric Binary Waveform. 

The restricted 2 PN waveform polarisations for non-spinning eccentric binaries are given by 
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The components (r(i), <I>(i)) denote the orbital seperation and phase of the system (also referred to as the true 
anomaly), and r{t) — dr/dt and $(i) = d^/dt. We noted that in the circular binary case, the angular frequency uj 
is defined hy uj ~ d^/dt. This is not true for eccentric binaries and furthermore, d^/dt is no longer a monotonic 
function of time. We will however continue to define the PN velocity parameter x as a; = (muj)'^'^ , but now where 
(jj = (27r + A<I>)/P = n + A$/P. Here n is the mean motion of the binary system and P = 2Tr/n is the radial orbital 
period. We should point out here that P is defined as the time to go from pericenter to pericenter. Due to precession 
effects, this is different from the time taken to go from $ to $ + 27r. The parameter A$ represents the advance of 
the pericenter per period. It was shown in Ref [58] that while the dominant spectral component for a circular binary 
appears at fgw = n/ir, this changes to (1 + k)n/-K for an eccentric binary having a PN accurate orbital motion, where 
k = A$/27r. 

Previous works evolve the eccentric binary in terms of the parameters (n, e) |22l. |58I462| | . where again n is the mean 
motion and e is the eccentricity of the system. However, Hinder et al [63| demonstrated that using the parameters 
{x, e) provides a better match with waveforms from numerical relativity. For this reason, we will also work with the 
parameter pair (a;, e). To go from {n, e) to {x, e) we use the following 2 PN relation 
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It was shown that eccentricity effects the waveform of a SMBHB in three ways [23[. First of all, the eccentricity 
induces amplitude modulation, secondly, it increases the amplitude of the waveform, and thirdly, it decreases the 
duration of the signal. 



The PN Conservative Orbital Dynamics. 



In order to write the relevant quan tities in terms of [x, e), we substitute the right hand side of the above equation 
into the necessary expressions in [61| which are expressed in terms of (tt., e) and truncate the subsequent equations at 
the 2 PN order. The parameters describing the conservative orbital dynamics can now be described in terms of the 
quantities {x{t),e{t),u{t)) 
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In the above equations, l{t) is the mean anomaly and is found by integrating the foUowing identity 
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The quantity u{t) is the eccentric anomaly and is determined by solving the transcendental 2 PN Kepler equation 
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We should point here that one of the differences in the evolution of the orbital phase between circular and eccentric 
waveforms, is that part of the secular evolution of $(t) for an eccentric binary is due to advance of periastron. This 
effect is taken into account by the leading term in the expression for $(i) given above. 



2. The PN Radiation Reaction 

To describe the effects of radiation reaction, the adiabatic evohition of x{t) and e(i) are given by the set of coupled 
1st order differential equations 
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These quantities provide the dominant secular evolution for the conservative orbital parameters (r, f, $, $, /). We can 
see that the above equations are functions of (a;, e) only. In order to generate the polarizations of the waveform, we 
first evolve this set of coupled ODEs. Once we have (a:(i), e(t)), we can integrate Eqn pO| for l{t), and then solve 
the transcendental Kepler Equation, i.e. Eqn (PT|) . for u{t). We then have all necessary quantities to then evolve the 
conservative orbital parameters. 

The two functions k_e(6) and Kj(e) which appear in the expressions for the tail terms at the 1.5 PN order in (i, e) 
are defined in terms of the following infinite series 
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where Jp{pe) are Bessel functions of the first kind, and JLipe) is given by 

J'{pe)^Jp-,{pe)^^-^^. 
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Rather than evaluate the summations themselves, we expanded both expressions in terms of truncated shifted Cheby- 
shev series 



K,E,j ~^akTl{y) 



(30) 



fc=0 



where Tj^{y) denotes the shifted Chebyshev polynomial over the domain 10 ^ < e < 0.6, and where 
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It was shown in |64| that the truncation error in a shifted Chebyshev series is equal to the coefficient of the truncated 
term. In this case, demanding a maximum error of 10"^" over the required interval, the infinite series for KE{e) and 
Kj(e) reduce to the following finite shifted Chebyshev series 
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We have also verified that outside of these intervals, the Cheybshev series have 
the shifted Chebyshev polynomials, we can use the fact that 

T*{y) = l, T*{y)^y, 

which now allow us to use the recurrence relation 
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to calculate the higher order shifted polynomials. 

We should point out that in its formal presentation, the above waveform description does not have support at 
e = 0. This can be seen from the fact that at e = 0, /S^ = oo and {u — v) is undefined. On the other hand, ke 
and Kj asymptotically approach unity as e — > 0, but again are undefined at e = due to infinities appearing in the 
expressions. We should also mention that at values of eccentricity of e < 10^^ the formal expressions for ke and kj 
become numerically unstable and begin to oscillate around unity. While our goal was always to compare the circular 
templates that have been used in the literature with eccentric templates, it is for these reasons that we need to use 
specifically circular templates and not just set e = in the eccentric templates. 



IV. ASTROPHYSICAL MODELING OF LOW AND HIGH MASS SEED BLACK HOLE BINARIES. 

To produce trustworthy catalogs of coalescing eccentric MBHBs, we proceed in two steps. Firstly, we extract 
catalogs of MBHB masses and redshifts from standard models of hierarchical MBH formation and evolution; then, 
we track the eccentricity of each system applying an hybrid model for the binary evolution in stellar dominated 
environments. 

A. Cosmological population of massive black holes 

MBHs are a ubiquitous components of nearby galaxy nuclei [see, e.g., |65|, and their masses tightly correlate with 
the properties of the host 1661 and reference therein]. In popular ACDM cosmologies, structure formation proceeds 
in a hierarchical fashion [67|, in which massive galaxies are the result of several merging events involving smaller 
building blocks. In this framework, the MBHs we see in galaxies today are expected to be the natural end-product of a 
complex evolutionary path, in which black holes seeded in proto-galaxies at high redshift grow along the cosmic history 
through a sequence of mergers and accretion episodes [48, 68] . Hierarchical models for MBH evolution, associating 
quasar activity to gas-fueled accretion following mergers between galaxies, have been successful in reproducing several 
properties of the observed Universe, such as the present day mass density of nuclear MBHs and the optical and X-ray 
luminosity functions of quasars [4^, |60] . 

In this general picture, the mechanism responsible for the formation of the first seed BHs is not well understood, 
and two distinctive families of models have became popular in the last decade. In the first family, seeds are light 
[M>100Mq, 'light seed' scenario |48|. being the remnant of the first POPIII star explosions [701; i^ ^^e second 
one, alre ady quite heavy (M>10^Mo) seed BHs form by direct collapse of massive proto-galactic discs ['heavy seed' 
scenario l49l. iTlj . The two models adopted here, representative of the two scenarios, are the same employed by the 
LISA parameter estimation task force 1131 (for the reader interested in more details, our light and heavy seed models 
correspond to models LE and SE in [T^l, respectively). Each model is constructed tracing backwards the merger 
hierarchy of 220 dark matter halos in the mass range 10^^ — 10^^ M© up to z = 20, using an extended Press & 
Schechter (EPS) algorithm (see [4^ for details). The halos are then seeded with seed black holes and their evolution 
is tracked to the present time. Following a major merger (defined as a merger between two halos with mass ratio 
M2/M1 > 0.1, being M2 the mass of the lighter halo) MBHs accrete efficiently, at the Eddington rate, an amount of 
mass that scales with the fifth power of the host halo circular velocity, normalized to reproduce the observed local 
correlation between MBH mass and the bulge stellar velocity dispersion (M — a relation, see [73 and references 
therein). For each of the 220 halos, all the coalescence events happening during the cosmic history are collected. The 
outputs are then weighted using the EPS halo mass function and integrated over the observable volume shell at every 
redshift to obtain numerically the coalescence rate of MBHBs as a function of black hole masses and redshift (see, 
e.g.. Fig. 1 in [zj)- In other words, the outcome of this procedure is the numerical distribution d^ N / dzdmidm2dt. 
We then perform 1000 Monte Carlo sampling of the d'^ N / dzdmidm2dt generated by each model, producing 1000 
catalogues of coalescing binaries over a period of one year. Distribution of MBH masses and mass ratios predicted by 
the two models are plotted in figure [T] The mass distribution is peaked around 10^ M0 for the light seed scenario, and 
around 10^ M© for the heavy seed scenario. Both mass ratio distributions are peaked at q — 1112 /rni = 1, however, 
the low mass seed model has a more gentle behavior, and coalescences are spread in the mass ratio range 0.1 — 1. 

B. Dynamics of eccentric massive black hole binaries in stellar environments 

The next step is to attach to each binary in the catalog an orbital eccentricity at some selected frequency, relevant 
for GW observations. Here is the description of our methodology. Following a major merger, the two MBHs sink to the 
center of the new formed galaxy because of dynamical friction. In star-dominated environment, N-body simulations 
pll |42| have shown that dynamical friction is efficient in driving the two MBHs down to a separation at which the 
enclosed stellar mass in the binary is of the order of TO2, without significantly affecting the stellar density profile. This 
is an indication that the evolution is still driven by the dynamical friction exerted by the overall distribution of stars, 
rather than by close individual encounters with stars intersecting the binary orbit. In our model, we assume that 
TO2 is driven by dynamical friction down to a separation uq, where the enclosed stellar mass in the binary is twice 
the mass of the secondary MBH. At that point, we apply the hybrid model developed by Sesana [43] to follow the 
semimajor axis and eccentricity evolution of the MBHB, assuming that the dynamics is purely driven by interactions 
with stars, and the stellar density proffie is characterized by a double power law pr = p{ri){r/ri)~'^ , with 7 = 2 for 
r > ri and 7 = 1.5 for r < r^. Here r^ = (3 — j)G{mi + m2)/cr^ is the influence radius of the MBHB, and p(r-,;) is the 
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FIG. 1: Mass (right panel) and mass ratio (left panel) distribution of coalescing MBHBs according to our MBH evolution 
models. The solid and the dashed histograms refer to the heavy and the light seed model respectively. The histograms are 
normaUzed so that the integrals in dlogq and dlogm are unity. 



stellar density normalization at r^. The reader is referred to [47] for full details, in the following we summarize the 
main features of the model. 

On its way to final coalescence starting from oo, the binary is subject to three main dynamical mechanisms driving 
its evolution. 

(i) Erosion of the cusp bound to the primary MBH. In this early stage, the MBHB extracts energy and angular 
momentum from the stars bound to the primary hole. During this process, lasting 10^ — lO^yr depending on the 
details of the system, the binary shrinks by a factor of ten, and e usually increases by a large factor, depending on 
the binary mass ratio and on the cusp slope. Results are tabulated in Table 1 of 40]. 

(ii) Scattering of unbound stars supplied into the binary loss cone by relaxation processes. After the cusp has been 
depopulated, further hardening is provided by super-elastic scattering of unbound stars diffused into the so called 
binary loss cone [36]. In general, hardenin g by scattering of unbound stars becomes effective when the binary reaches 
the so called hardening radius, defined as [33] ah ~ Gr?i2/(4tT^). This is the separation at which the specific binding 
energy of the binary is of the order of the specific kinetic energy of the field stars. Once the binary is hard, its 
hardening proceeds at about constant rate and,^hile circular binary tends to stay circular, even slightly eccentric 
binaries tend to increase their eccentricity f37[ |39| . Whether the binary reaches the point at which GW emission 
becomes efficient (typically a gw5i lO~^Qft., see [43]), depends on the rate at which stars are supplied to the MBHB loss 
cone. However, as noted by [47] the typical eccentricity evolution of the system depends only mildly on such rate. 
Here we assume efficiently repopulated (also referred to as 'full') loss cone. In this case the binary typically enter the 
GW-dominated phase in a timescale of <10*yr. 

(iii) Emission of GWs. The effect of GW emission is modeled in the quadrupole approximation up to a selected 
typical GW frequency (/gw = 2/^, being fk the Keplerian frequency of the binary), which is 2 and 6 x 10~^ Hz for 
the high and low seed sources respectively, corresponding to orbital separations of many lO's to lOO's of m. Under 
this assumption, the evolution equations for the system are given by [20] 
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FIG. 2: MBHB evolutionary tracks produced by our model for selected MBHBs. In the top and central panels we plot the 
evolution of e and a/ao as a function of time, respectively; in the lower panels we represent the evolution of e as a function of 
the Keplerian frequency fk of the system. Different line-styles refer to eo = (solid), 0.3 (short-dashed), and 0.6 (long-dashed). 
Dotted vertical lines in the lower panels mark fk = few/2 — 10~^ and 3 x 10~^ Hz. 

The shrinking rate is a strong factor of a, meaning that GW-driven hardening is effective only at small separations. 
The eccentricity evolution rate is also a strong function of a and e itself, and it is always negative. GW emission, 
therefore, is very effective in circularizing MBHBs, which, in turn, is the reason why little attention has been paid so 
far to eccentric systems in the context of GW detection. 

We numerically solve two coupled differential equations for the evolution of a and e, by combining the three 
dynamical mechanisms mentioned above, as detailed in Section 2.3 of [4^. We consider three different values of 
Co = 0,0.3,0.6, at the initial semimajor axis ap. This accounts for the fact that galaxies typically capture each other 
on very eccentric orbits, which are reflected in the initial trajectory of the two MBHs. Even though dynamical friction 
against massive gaseous disks has been proven efficient in circularizing the orbits of the two MBHs [34], this is in 
general not true in stellar dominated environments 26], and the two MBHs may pair together on a significantly 
eccentric orbit at ap. 

In figure^ we plot the evolution of the MBHB semi- major axis and eccentricity, for values of ttii and q representative 
of both the light and the heavy seed scenarios, for different values of eg. As discussed in [43], when cq — 0, the 
eccentricity growth is more significant for unequal mass MBHBs ; while binaries with eo>0.3 reach values of e close 
to unity irrespective on the other binary parameters. Lighter binaries generally preserve higher residual eccentricities 
in the LISA band. The lower panels represent the evolution of e versus the orbital Keplerian frequency fk- When we 
start the PN evolution (fk = /gw/2 = 10~^ and 3 x 10^'"' Hz for the heavy and the light seed models, respectively), 
MBHB eccentricities are in general less than 0.1 for MBHBs with ep = 0, but they can be as high as ^ 0.8 when 
eo>0.3. 

We now have six catalogs of MBHBs, three for each seed model, assuming eo = 0, 0.3, 0.6. Each catalog contains a 
Monte Carlo realization of the coalescing MBHB population and includes the source redshift z, mi, m2, and e at the 
selected reference Keplerian frequency. 



COMPARING CIRCULAR AND ECCENTRIC TEMPLATES FOR LISA. 



As previously stated, our aim in this study is to evaluate the fidelity of restricted PN circular inspiral templates in 
the presence of eccentric SMBHBs. A full exploration of eccentric black hole binary parameter estimation is currently 
underway 75|, but is outside the scope of this particular study. To compare the template families, we ran a Monte 
Carlo simulation where we assumed a one year LISA mission lifetime. The main reason for this is, most of the SNR 
comes from the end of the inspiral. Therefore, as the starting point of observation is arbitrary, it is our ability to 
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FIG. 3: The source frame individual mass distributions as a function of redshift for low mass (top) and high mass (bottom) 
seed black hole binaries. The three cells going from left to right, both top and bottom, show the catalog individual masses 
assuming initial eccentricities at binary hardening of eo = 0, 0.3 and 0.6 respectively. In all cells m\ axe denoted by the black 
circles, and m2 are represented by the red squares. 

match the eccentric waveform for a coalescing binary in the final year that is important. We also investigate the low 
and high mass seed cases individually. We will explain the organization of the Monte Carlo runs in greater detail 
below. 



A. The Monte Carlo Setup. 



The first point of importance is termination of the eccentric inspiral waveforms. As we do not know the exact 
position of the last stable orbit (LSO) for eccentric binaries, we decided to terminate the waveforms before the orbital 
separation reached r — 6m (i.e. the LSO for a test particle in a Schwarzschild geometry) or 5 mHz, whichever is 
reached first. In Figures ^ and ([5]), we plot the evolution of eccentricity at -3 years, -1 year and at the final point 
of evolution. These catalogues provide the initial conditions in terms of eccentricity for our simulations. 

The Monte Carlo was set up as follows : We started with six source catalogues consisting of high or low mass seeds 
with initial eccentricities of eo = 0, 0.3 or 0.6 at binary hardening. For each of the sources in the eccentric catalogues, 
we randomized the sources over the parameters (t, ip, ipc, 0, 0), while using the catalogue values of (toi, ?ti2, z, xq, eo). 
For each system, we calculated the optimal signal to noise ratio (SNR) to ensure we achieve a threshold value of 8. 
This value was chosen as it is a value at which many algorithms start to discern between a real SMBHB signal and a 
signal produced by the multitude of white dwarf binaries in the galactic foreground. The optimal SNR is defined as 



where h,, denotes an eccentric template, and the angular brackets represent the noise weighted inner product 

{a\b) =2 r ^~a[f)b*{f)+cc, 



where 



«(/) 






(38) 



(39) 



(40) 



is the Fourier transform of the time domain function a{t) and the cc denotes complex conjugate. The integration 
limits {fQ,fh) are given by /o = 10^^ Hz, and fh — fgwir — 6m) or 5 x lO^'^ Hz, depending on the source. The 
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FIG. 4: Eccentricity evolution over the final three year period for low mass seed black hole binaries. The plot displays the 
evolution for the three models of initial eccentricity. Note that the evolution reads downwards for each model. 
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FIG. 5: Eccentricity evolution over the final three year period for high mass seed black hole binaries. The plot displays the 
evolution for the three models of initial eccentricity. Note that the evolution reads downwards for each model. 



quantity 5'„(/) is the one-sided power spectral density and is composed of a combination of both instrumental noise 
and the confusion limited galactic background, i.e. 



Sn{f)^ST''{f)+Sr'{f). 



where the instrumental noise S^^^^^{f) is given by 
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and the confusion noise estimate Sf° (/) is derived from a Nelemans, Yungelson, Zwart galactic foreground model ( J76l . 
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In the above expressions, L = 5 x 10'' km is the arm-length for LISA, S'p°*(/) = 4 x 10"^^m2/i72 and S^'^'^f) = 
9 X 10~30;^2^g4^^^ g^j.g iijg position and acceleration noises respectively. The quantity /* — 1/{2ttL) is the mean 
transfer frequency for the LISA arm. The instrumental noise also contains a reddened noise term which steepens the 
noise curve between lO""*^ and 10~^ Hz. Finally, the units of 5'™"^(/) are WxT^ . 

If the optimal SNR is greater than the threshold of 8, we then run a short Markov Chain Monte Carlo (MCMC) 
with circular templates, using the exact parameters of the eccentric binary system as our initial best guess for the 
circular templates. The MCMC is a stochastic search method which has been used many times in LISA GW data 
analysis. While these works used sophisticated variants of the MCMC to search over a wide parameter space, in this 
study we arc only concerned with local exploration and are therefore using a straightforward MCMC. While we refer 
the reader to 2] for an indcpth discussion about MCMC methodology in GW data analysis, briefly, the method works 
as follows : starting with the signal s{t) = he{t) +n{t), where n(t) is the noise in the detector, and some initial circular 
template hc(t; Xi) constructed by choosing a random starting point in the parameter space Ai, we then draw from a 
proposal distribution and propose a jump to another point in the space A,;+i. In order to compare both points, we 
evaluate the Metropolis-Hastings ratio 



H = 



n{Xi+i)p{s\Xi+i)q{Xi\Xi+i) 



n{Xi)p{s\Xi)q{Xt+i\Xi) 
Here 7r(Ai) are the priors of the parameters and p{s\Xi) is the likelihood defined by 

£(a,) =Ce-(^-"=(^0N-''c(^0>/2, 



(43) 



(44) 



where C is a normalization constant. The quantity q{Xi\y) is the proposal distribution used for jumping from Xi to 
Ai+i. For this study, the proposal distribution is a multivariate Gaussian calculated using the Fisher information 
matrix (FIM) 



/ dhg 



dX" 



(45) 



Once a jump is proposed, it is then accepted with probability a = min{l, H), otherwise the chain stays at A^. 

In order to speed up the convergence of a Markov chain, it has been shown ^ that heating the likelihood surface 
via simulated annealing helps the chain to move more easily by effectively smoothcning and reducing the height of 
maxima on the surface. To this end, we run the first 1000 iterations of the MCMC using a simulated annealing phase 
which replaces the value of 1/2 in the exponent of Eqn (pi)) by a factor /3 where 



ilO-«(i-*) 0<i<Tr 



(46) 



i>Tr, 



For this particular study, ^ is the heat-index defining the initial heat and is taken to be equivalent to the optimal SNR 
calculated using the circular templates, i is the number of steps in the chain and Tc = 1000 is the cooling schedule. 
We should remark here on our choice of the initial heat. We tried scaling the heat to the optimal SNR for eccentric 
binaries, but this make the initial heat very high and as we start the MCMC close to the true solution, this amounts 
to a wasted number of computer cycles. As in general, the optimal SNR of a circular binary is less than the optimal 
SNR of an eccentric binary (something we will justify at a later point), we use the circular binary SNR to scale the 
initial heat. 
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One of the quantities that we investigate is the overlap between the template and the signal. While the concept of 
a global overlap for LISA is ill-defined, the overlap in a particular channel is defined by 

= —^d^^=. (47) 

^{hc I he) {he I he) 

While the MCMC changes the parameters of the system, the overlap is a simple way of measuring the improvement 
in the fit between the eccentric and circular templates. 

VI. RESULTS. 

In the following sections we present the results for the 1 year missions for both the low and high mass seeds 
seperately. 

A. High Mass Seed Black Hole Binaries 

In Fig ([SJ we plot both the optimal and maximum recovered SNRs from the high mass seed Monte Carlo's. We can 
see from the top panels that optimal SNRs for eccentric binaries peak at values of 50-100, with varying maximal SNRs 
of between 600 and 800. In the bottom panels we have plotted the maximum SNR recovered by a maximized circular 
template at the end of the MCMC. We can see that for all eccentric catalogues, the distributions peak at SNRs of 
5-10, and have maximum values of 25-30. This means that for the high mass seed systems, the circular templates 
have median optimal SNR recovery of just 5, 3 and 5% for the three models of initial eccentricity. 

To properly explain this result, we refer to Fig ([7]). In the top panels we plot the initial overlap between a circular 
and eccentric waveform with exactly the same parameters in both the A (blue) and E (red) LISA channels. We can 
see that while the distributions vary from slightly negative to slightly positive, they are roughly peaked around zero. 
This implies that circular and eccentric templates with exactly the same parameters are essentially orthogonal to each 
other, regardless of the value of the initial eccentricity. In the bottom panels we show the distributions of maximum 
overlaps at the end of the MCMC. It is clear that the MCMC has managed in each case to find a more suitable 
parameter set in terms of improving the likelihood between the circular and eccentric templates, but the overlaps are 
no-where near what would be required for LISA data analysis. The best that we could do with the circular templates 
is confirm a detection, as long as the threshold for a detection is sufficiently high. 

To further explain the disparity between circular and eccentric templates, we have plotted the spectra of a some 
SMBHB systems in Fig ([5]). The system in the left hand panel has individual redshifted masses of 4.238 x 10^ and 
3.728 X 10^ Mq at z = 4.49. For this system the initial and final eccentricities are e^ = 1.8 x 10"'* and e/ — 2.3 x 10~^. 
The system in the right hand panel has individual redshifted masses of 1.975 x 10^ and 2.157 x 10^ M0 at z = 6.9. 
For this system the initial (this we define as being the eccentricity at the beginning of the observation) and final 
(when the binary separation reaches 6m or 5 mHz) eccentricities are e^ = 0.438 and e/ = 0.011. Both systems have 
the same angular parameters and sky locations. If we first focus on the system in the left hand cell, while it looks 
like there is a good match between the eccentric (blue) and circular (red) templates at low frequencies (we should 
point out here that if we zoom in on these frequencies, there is a clear phase mismatch between the waveforms), it 
is clear that even for an extremely mildly eccentric binary, there is extra power at higher frequencies that we do not 
see in circular templates. As quantities such as overlaps and SNRs require good phase matching, we can now see why 
circular templates do such a bad job in capturing eccentric binaries for LISA. In the panel on the right, for an eccentric 
waveform with substantial eccentricity, it is clear that there is very little possibility of a circular template being able 
to match the higher power and frequency content of such an eccentric binary. Even if an algorithm managed to find a 
system with low enough masses that it pushed the LSO frequency close to that of the eccentric binary, it would still 
not be able to capture the structural information. 

While it will not be possible to carry out a parameter recovery of eccentric binaries using circular templates, as 
we said earlier, it may still be possible to use the circular templates in a detection only framework. Thus, it is 
still interesting to investigate the parameter mismatch between the two waveform families. In Fig. © we plot the 
parameter mismatch at the end of the MCMC for chirp-mass, reduced mass, luminosity distance and sky position. 
For the sky position error Atr, we calculate the orthodromic distance between the sky position of the eccentric binary 
and the sky position of the best fit circular template. This is done using a special case of the Vincenty formula (which 
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FIG. 6: The top row of this figure shows the optimal SNR distribution using high mass seed eccentric waveforms for each of 
the three initial eccentricities. The bottom row shows the maximum SNRs recovered by the MCMC using circular templates. 
We can see a massive degradation in the recovered SNR using circular templates. 
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FIG. 7: The top row of this figure shows the initial distribution of overlaps achieved by comparing circular and eccentric 
templates for high mass seeds with exactly the same parameters for each of the three initial eccentricities in the LISA A 
(blue) and E (red) channels. The bottom row shows the distribution of overlaps at the end of the MCMC. While there is an 
improvement in fit, the overlaps are not close to what we require for LISA. 

is normally used for calculating the distance between two points on an ellipsoid): 



Act = arctan 



(cos (pr sin A6l) + (cos (f>a sin (j)r — sin (pa cos (j)r cos A9l) 

sin (pa sin (p^ + cos (pa cos (pr COS A6l 



(48) 



where Act has units of radians, {(paj^f'r) denote the actual and recovered longitudes and A6l = 01 ^ ^l i^ ^^^ 
difference in actual and recovered latitudes. We have used this particular expression as to avoid the large rounding 
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FIG. 8: A comparison of the spectra of two circular and eccentric systems with identical paramters. The eccentric system on 
the left corresponds to a high mass seed system with initial and final eccentricities of 1.8 x 10~* and 2.3 x 10~^, while the 
system on the right has initial and final eccentricities of 0.438 and 0.011 . We can see that in both cases the circular templates 
are not a good match to the eccentric waveforms. 

errors associated with the spherical law of cosines for the case where Ac <c 1 and with the haversine formula in the 
case of antipodal points. 

We can see from Fig ([9]) that for all three initial eccentricity models, the errors in the parameter estimation are 
bigger than we are used to for these types of sources [2]. If we first focus on the mass parameters, for the models 
where the initial eccentricity was or 0.3, we can see that the chirp and reduced mass fractional error distributions 
are peaked at approximately 10~^ and 10~^ respectively. As a comparison, MCMC searches have recovered fractional 
errors in both mass parameters on the orders of 10~^ and 10~^ for circular SMBHBs 2]. For the cq = 0.6 case, 
we can see that while a number of the sources are resolvable with similar precision as in the other two cases, there 
are a number of sources where the mass parameters are essentially undetermined. Now focusing on the error in the 
sky position, for all three cases, we can see that the sky is essentially undetermined with massive errors in the final 
estimated sky position. Finally, for the estimation of luminosity distance, we end up in a situation where most of 
the distances are unresolved. The overwhelming conclusion here is, while for the higher SNR cases we may have 
confidence in a detection, we can not have confidence in the system parameters extracted using circular templates. 



B. Low Mass Seed Black Hole Binaries. 



For low mass seed SMBHBs, the situation is slightly better. There are a few reasons for this : as we previously 
stated, if a system does not reach the orbital separation of r = Qm before reaching a GW frequency of 5 mHz, we 
terminate the waveform at 5 mHz. Because of the masses involved, the coalescence frequencies of many of the systems 
are at lO's of mHz and are thus outside of our band of interest for this study. This effect can be seen in the optimal 
SNRs presented in the top panels of Fig (ITUl) . While we again have SNRs of 400-600, the distributions are peaked at 
SNRs of 10-20. We know from previous studies that most of the SNR is recovered from the final few cycles, usually 
corresponding to the last few days of inspiral. As we do not see the coalescences for most of these systems, the circular 
templates have to fit a less relativistic waveform and thus have the ability to perform better. Also, we can see from 
Figs (HJ and ([5]) that the eccentricity evolution is slower for the low mass seeds. This means that while the eccentricity 
of low mass seeds stays higher for a greater period of time, there is more possibility for the circular template to obtain 
a better fit as the system is not changing as quickly as the high mass seed case. Finally, we remarked earlier that 
eccentric binaries radiate at higher harmonics of the orbital phase. In the low mass seed case, as we are truncated 
virtually all of the systems at 5 mHz, the effect of the radiation at higher harmonics is not as influential as in the 
high mass seed case (we do note, of course, that this is due to the use mainly to our use of the LFA for the LISA 
response). This effect has also been seen in the case of circular binaries where higher harmonic corrections have been 
added [il[.We can see from the bottom panels of Fig P^ that while the circular templates still suffer, they actually 
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FIG. 9: High mass seed black hole parameter mismatch for chirp-mass, reduced mass, sky position and luminosity distance for 
the three models with initial eccentricity of eo = (solid line), eo = 0.3 (dashed line) and eo = 0.6 (dot-dashed line). 

recover a greater percentage of the optimal SNR than in the high mass seed case. While the SNR distributions are 
still peaking at low values of 5-10, we are seeing systems with maximum recovered SNRs of between 40 and 80. For 
the low mass seeds, we are now achieving median optimal SNR recoveries of 13, 10 and 12% respectively for the three 
models of initial eccentricity. 

If we now move our attention to Fig (jlip. we again plot the initial overlaps for circular and eccentric templates 
with identical parameters for the LISA A (blue) and E (red) channels in the top panels. We see a similar story here 
to the high mass seed case where the initial overlaps are peaked around zero, again showing that the two template 
families with identical parameters are essentially orthogonal to each other. In the bottom panels we plot the maximum 
overlaps at the end of the MCMC. Again, while the peak of the distributions is still close to zero, we do see systems 
with overlaps approaching 0.25, showing an improvement over the high mass seed case. We should once more point 
out that with identical parameters, the initial overlaps seem to be independent of the initial eccentricity. 

While a maximum overlap of 0.25 is an improvement, it would not really give us enough confidence in our detection. 
However, just as in the high mass seed case, it is interesting to look at the effect of parameter estimation in the low 
mass seed case. In Fig (I12p we plot the errors in parameter estimation for the chirp and reduced masses, the sky 
position and the luminosity distance. It is here that we see the benefit of not having to fit the merger of the waveform. 
We can now see that for all three models, the fractional errors in the chirp mass estimate peak somewhere between 
10~^ and 10~^. For the reduced mass, while some systems are unresolvable in this parameter, the vast majority peak 
with errors between 10~^ and 10~^. However, this is where the good news ends. We again see that in all three models 
both the sky position and luminosity distance are essentially unresolved for all systems. So once again, the conclusion 
is that also in the case of low mass seeds, circular templates are not efficient enough in capturing eccentric black hole 
binaries. 



VII. CONCLUSION. 



In this work, we have used a number of source catalogues taken from the end stage of a hybrid model for the 
evolution of eccentric SMBHBs. These catalogues describe both high and low mass seed systems with initial model 
eccentricities of 0, 0.3 and 0.6. These systems were then evolved into the LISA detection window using the Peters and 
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FIG. 10: The top row of this figure shows the optimal SNR distribution using low mass seed eccentric waveforms for each of 
the three initial eccentricities. The bottom row shows the maximum SNRs recovered by the MCMC using circular templates. 
We can see a massive degradation in the recovered SNR using circular templates. 
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FIG. 11: The top row of this figure shows the initial distribution of overlaps achieved by comparing circular and eccentric 
templates for low mass seeds with exactly the same parameters for each of the three initial eccentricities in the LISA A 
(blue) and E (red) channels. The bottom row shows the distribution of overlaps at the end of the MCMC. While there is an 
improvement in fit, the overlaps are not close to what we require for LISA. 



Mathews' equations to model the secular decay of the semi-major axis a and the eccentricity e due to GW emission. 
At a certain point, we then changed over to the 2-PN equations for the evolution of the PN velocity parameter x and 
eccentricity. This evolution then provided us with the initial conditions to examine the sources in the final year of 
evolution, either to coalescence or a maximum GW frequency of 5 mHz. 

An extensive study was then carried out by combining a Monte Carlo simulation over the extrinsic parameters of 
the system, combined with a Markov chain Monte Carlo, to examine the fidelity of searching for eccentric systems 
using circular templates. We found that in the high mass seed case, only about 5% of the optimal SNR was recovered 
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FIG. 12: Low mass seed parameter mismatch for chirp-mass, reduced mass, sky position and luminosity distance for the three 
models with initial eccentricity of eo = (solid line), eo = 0.3 (dashed line) and eo = 0.6 (dot-dashed line). 



by the circular templates. Worse still was that the fit between the eccentric and circular waveforms at the end of 
the MCMC only achieved overlaps of about 0.15, much below the confidence level needed for LISA data analysis. 
We also looked at the errors in parameter estimation and found that while the mass parameters were resolvable, the 
errors were quite large. For all three initial eccentricity models, both the sky position and luminosity distance were 
unresolvable. This has the consequence that as the redshifted masses are similar to the true redshifted masses, using 
a particular cosmological model, we would interpret the true system as being composed of much lighter binaries at a 
higher redshift. 

For the low mass systems, the recovered optimal SNR was better with a maximum of 13%. But again, the maximum 
overlaps were only on the order of 0.25. While this represented an improvement over the high mass seed case, and 
can be attributed to the fact that for the vast majority of these systems we did not see the coalescence, thus meaning 
that the circular templates had to fit a less relativistic waveform, it is still far below the LISA confidence level. While 
the errors in the estimation of both the chirp and reduced masses were smaller in this case, the sky position and 
luminosity distance were again unresolved. 

The clear consequence of this study is that LISA data analysis will require the use of eccentric templates, even if the 
final eccentricity is on the order of 10^^ < e/ < 10""*. While ground based studies have demonstrated that circular 
templates work at this level, the case for LISA is very different. It was found for LIGO that fitting factors were 
reduced when neutron star binaries were considered, rather than black hole binaries. This was due to the fact that a 
stellar mass black hole binary lasts less than one second in the detector, whereas a neutron star binary is observable 
for a little over 20 seconds. In the LISA case, with signals lasting many months or years, there is much more scope 
for the circular and eccentric templates to be out of phase with each other due to the total number of cycles, so our 
overall result is really not surprising. 

Our study also has consequences for Numerical Relativity as a lot of work has gone into reducing the residual 
eccentricity in merger waveforms. If it is true that eccentric binaries will be observable in LISA, with some of the high 
eccentricities seen in this study, then we will require a catalogue of merger waveforms for eccentric binaries before 
LISA launches. 

We have now begun two further studies regarding eccentric binaries : the first is an extensive parameter estimation 
study. In the second study, we intend to revisit the current problem but using circular binaries with higher harmonic 
corrections. 
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